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Abstract. - The total energy E(t) in a fluid of inelastic particles is dissipated through inelastic 
collisions. When such systems are prepared in a homogeneous initial state and evolve undriven, 
E(t) decays initially as t~ 2 ~ exp[— 2er] (known as Haff's law), where r is the average number of 
collisions suffered by a particle within time t, and e = 1 — a 2 measures the degree of inelasticity, 
with a the coefficient of normal restitution. This decay law is extended for large times to E(t) ~ 
r~ d / 2 in d-dimensions, far into the nonlinear clustering regime. The theoretical predictions are 
quantitatively confirmed by computer simulations, and holds for small to moderate inelasticities 
with 0.6 < a < 1. 



A fluid of inelastic hard spheres (IHS) in d dimensions, prepared initially in a state of 
thermal equilibrium with temperature To, will remain for an extended period of time (measured 
in average number of collisions suffered per particle) in a spatially homogeneous cooling state 
(HCS), provided the degree of inelasticity e is small. In this state the average kinetic energy 
per particle E(t) = %T(t) is decreasing like 1/t 2 due to inelastic collisions, as first explained by 
Haff [Q. However, the HCS with a spatially uniform density and temperature, and a vanishing 
flow field is unstable against long wavelength spatial fluctuations, and leads finally to a state 
with large scale clusters in the density field n(r, t) as well as in the flow field u(r, t) (vortices). 
This is illustrated in the MD simulations of fig. [IJ In this state cooling slows down and 
large deviations from Haff's law occur, which have not yet been explained in any quantitative 
manner j| |j. 

In general analytic results for nonlinear decay of rapid granular flows are rare The 
goal of this letter is to calculate the long time decay of E(t) from nonlinear hydrodynamics, 
using mode coupling methods. The basic idea is to decompose the quantities of interest, here 
the total microscopic energy, or — in the study of long time tails of Green-Kubo formulas — the 
total microscopic momentum-, energy- or particle fluxes, into a superposition of products of 
slow hydrodynamic fluctuations of density n(r, t), temperature T(r, t) and flow field u(r, t). 
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Fig. 1. - 

Snapshots (a) of flow field u(r, r = 80) and (b) of density field rt(r,r = f60) in an 1HS 
simulation with N = 50000, a = 0.9, <f> = 0.245, prepared in an initial thermal equilibrium 
state. Crossover to the nonlinear regime occurs at r c ~ 70. 

As one is dealing with fluctuations their amplitude is small in general, so their time decay can 
be calculated from linear hydrodynamic modes Ref. [Q, provided the fluctuations are stable. 
However, in the undriven 1HS fluid the density fluctuations are growing exponentially like 
exp(er) (clustering instability) The redeeming feature in the resolution of this paradox is 
that the energy only depends on products of modes, in which exponentially growing density 
fluctuations are balanced by exponentially decreasing temperature fluctuations. 

In the case of inelastic hard spheres two colliding particles loose a fraction, e = 1 — a 2 , 
of their relative kinetic energy, where e is the degree of inelasticity and a the coefficient of 
normal restitution. Particles move more and more parallel after every collision, and their 
random motion (temperature) decreases; i.e. the spatial fluctuations at shorter wavelength 
(which are rapidly damped out through viscosity and heat conduction) are not replenished 
with energy from the more microscopic degrees of freedom, as in the elastic case, where 
fluctuations at all wavelengths are kept at thermal noise level, fed by the randomizing elastic 
collisions. So, the dynamics selectively suppresses the shorter wave length components of the 
flow field, as is clearly illustrated in fig. la. After a time on the order of the mean free time 
to between collisions, the surviving fraction of the energy E(t) = (1/iV) (XZi \ mv 1) * s totally 
stored in hydrodynamic modes (with wavenumbers k < 1, where k is the measured in units of 
inverse mean free path Iq) of the kinetic energy of the flow and the internal energy (granular 
temperature), 

E(t)^l J dr|i<p(r,i)u 2 (r,i)) + ^<n(r,i)T(r,t))|. (1) 

The subsequent decay of E(t) is controlled by two different k-regimes ||: (i) the elastic 
regime (e < k < 1) at short times (homogeneous cooling regime) with to < t < to/t, and (ii) 
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the dissipative regime (fc < e) at large times (t 3> to/e). 

In the first regime the decay is controlled by the HCS, where n(r, t) = n, u(r, t) = and a 
homogeneous temperature, decaying according to Haff 's cooling law, 

T ^ = n I T \l+ \2 = T ° cx P(- 2 7or), (2) 

[1 + lot /tor 

where Eq = ^Tq is the initial (equilibrium) energy per particle, 70 = e/2d, and to — l/u(To) is 
the mean free time with collision frequency w(T ) in the initial state. The second equality in 
eq. (||) defines the average number of collisions per particle within a time t. It is obtained by 
integrating dr = uj{T(t))dt, where ui(T) ~ \/T/lo with a temperature independent mean free 
path lo, given by the Enskog theory for a dense system of hard spheres (d = 2, 3). Haff 's law is 
illustrated in figs. 2a, b by the straight lines, labeled 1. The HCS is essentially an adiabatically 
changing equilibrium state, parametrized by {n, u = 0,T(t)}. 

To evaluate the contribution of the long wavelength fluctuations, we linearize the fields in 
cq. (d) around the HCS, where u(r,i), 5n(r,t) = n(r,t) - n and ST(r,t) = T(r,t) - T{t), 
are small (at least their combined effects; see redeeming features below), and consider Fourier 
modes. Retaining up to quadratic terms yields, 

E(t) = ^T(t) + ^ N J ^-^{ p {\u(k,r)\ 2 )+d{Sn(k,T)ST(-k,r)}, (3) 

where T(t) is given by Half's law. The first term in the integrand is proportional to the 
structure factors [(d — l)S±(k,r) + S\\(k,r)], and the second one to S n T(k,r), where the 
subscript || refers to the longitudinal component of u and the subscript _L to the transverse 
ones. 

According to the theory of ref. [3a, b], the fluctuations around the HCS are described by 
fluctuating hydrodynamic Langevin equations, and obey an 'adiabatic' fluctuation dissipation 
theorem. The results from the structure factors, Sij(k, r), obtained by numerically integrat- 
ing the fluctuating hydrodynamic equations of ref. [3b], could be integrated over the full 
hydrodynamic range of wave numbers to obtain the decay of E(t) for all times t > to- 

Here, however, we present an analytic method that enables us to calculate the asymptotic 
decay for t 3> to/e explicitly. For such times only modes in the dissipative regime with k <C e 
will contribute to the integrands in (||). In this regime, the dispersion relations for the modes, 
relative to the HCS, are well known [5,3c]. 

The structure factors for the (d — 1) transverse components u± a and the longitudinal one 
ui in the range k <§; 70 ~ e can be deduced from ref. [3b], and yield, 

V- l (\u ±a (k, r)| 2 ) = (T /p)e- 2A ^ (l + A ± fc 2 /7o) , (4) 

with a diffusivity = v/I^lo, where v = rj/p is the kinematic viscosity. The longitudinal 
structure factor is given by a similar expression with replaced by A11. The explicit form 
of the longitudinal diffusivity (given in ref. [3b] as A|| = 7o£ 2 Ao) i s n °t needed here. One can 
verify from the explicit expressions that A11 > Aj_, and even A|| 3> A^ for 70 — > 0. 

It should be noted that the dimensionless diffusivities, A^ and Ay, do not depend on the 
local temperature, since transport coefficients and collision frequencies are proportional to \pT. 

By inserting eq. (^) into (^) and performing the k-integrals, one easily finds that the 
contributions E uu (t) of the term (|w| 2 ) to eq. ||) decay as, 

T n f d-l 1 1 / 1 ^ d/2 



Euu{t) ~ 2< 1 (A ± )<V 2 + (A||) d / 2 J [sttt) ■ 
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Next, we consider the contributions E n T{t) in (||) from the n-T fluctuations. Calculation 
of the time dependence of (5n(k, r) and <5T(k, t) in the dissipative regime is in principle 
straightforward, but technically much more involved. One solves the initial value problem for 
the coupled linearized hydrodynamic equations by determining its eigenvalues (to C(fc 2 )-terms 
included), and eigenvectors (to O(fe)-terms). The result for k <§; e is, 

<5n(k,r) ~ — Cie^°- A i!' £2 ) T u / (k,0) + --- 
7o 

6T(k,r) ~ l Ac 2 e-^ +A ^ T Ul (k,0) + --- : (6) 
7o 

where the dots represent subdominant decaying terms, and C n (n = 1,2) are some constants. 
The technical details are available through ref . || . 

The first equation in (^J) exhibits the clustering instability. It shows that the density 
fluctuation is growing at an exponential rate. We are not aware of any publication on 
a freely evolving rapid granular flows or IHS fluids, that explicitly demonstrate how the 
clustering instability is driven through the coupling to the long wavelength longitudinal velocity 
fluctuations — except for the closely related nonlinear slaving mechanism between density and 
flow field, discussed by Goldhirsch et al. || — , whereas all other fluctuations remain bounded 
by their values in the initial thermal equilibrium state, i.e. the fluctuations in the flow field 
decay diffusively at rates A±k 2 and Ayfc 2 , and those in the temperature decay rapidly at a 
rate (70 + A||fc 2 ). 

The above discussion also implies that an incompressible granular flow (V • u = 0) would 
not exhibit a clustering instability, at least not in a linear stability analysis. Moreover, an 
initial density fluctuation <5n(k, 0) would not cause any clustering instability. Furthermore, 
we observe that (Sn(k, r)5T(— k, t)) remains bounded, because the exponential increase of the 
first factor is balanced by an exponential decrease of the second. The resulting integral yields 
E nT (t) cx ©(r-^ 2 " 1 ), which is a subleading asymptotic term when compared to (j|). 

It can be shown similarly that the three mode contribution E nuu (t), neglected in the 
transition from eq. ([l]) to (|j|) remains bounded, and decays as T~ rf , again a subleading 
correction. The dominant large time contribution E uu {t) is plotted in figs. 2a,b, where the 
lines labeled 2 and 3 represent the contributions of the transverse and longitudinal modes, 
respectively. The former contribution is a factor e 3 ~ 20 larger than the latter. 

The result (^) also shows that there is no balancing of unstable density fluctuations in the 
structure factor S nn (k,T) ~ (\Sn(k, t)\ 2 ) . Consequently, the predictions of the hydrodynamic 
Langevin equations for S nn breakdown at the crossover from the homogeneous cooling regime, 
as explained in ref. [3b], whereas those for S±, S\\ and S n T, used in this paper are expected 
to remain valid until far into the nonlinear clustering regime. 

In summary, in the decay of the total energy E(t), one can distinguish two different 
behaviors: homogeneous cooling for r < r c , where Bhcs(^) — 1^0 ex P( — %yo r ) is given by 
Haft's law, and a diffusive decay for r > r c , where E(t) ~ E uu (t) is solely determined 
by the flow field, all long wavelength components of which decay diffusively. The equality 
E uu = Encs defines the crossover time r c , typical values of which are listed in the figure 
captions and table 1. 

Before comparing the preceding results with MD simulations, a final comment about the 
relation between the 'internal' kinetic time r and the 'external' time t is needed: in the 
homogeneous cooling state the relation between both times is correctly predicted by r = 
(I/70) m (l + lot/to) m eq. (^). In the nonlinear clustering regime the relation between both 
times is not understood. The value of r(t), as measured in the computer simulations for r > r c , 
increases faster than those given by (2). For the cases corresponding to figs. fk,b the time r 
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becomes roughly linear in t for t larger than 80. 

In fig. ||a the theoretical prediction for E(t) is compared with computer simulations in a 
dense system (<fi = 0.245) of N = 50000 IHS particles at low dissipation (a = 0.9). The 
density or coverage is defined as <fi — jnNa 2 /L 2 , where a is the diameter of a disk. The 
agreement remains good until far into the nonlinear clustering regime with r ^> r c ~ 70 (see 
fig. |l|a,b). The density field n(r,i) contains small (r = 80) and large (t — 160; see fig. 0b) 
spatial inhomogeneities. The flow field shows both at r = 80 (see fig. [j]a) and at r = 160 
well-developed vortex patterns. To describe the crossover regime around r c one would have to 
evaluate E{t) in the total hydrodynamic time regime (t > to), by performing the numerical 
calculations described in the paragraph below eq. (||). The same good agreement is found at 
small inelasticities (a = 0.975,0.9) for all densities (<fi = 0.05,0.11,0.245,0.4) in the larger 
systems with N = 50000 particles, as summarized in Table 1. When the number of particles 
is reduced to N = 20000 and N = 5000, the deviations steadily increase, suggesting finite size 
effects. 




Fig. 2. - 

Energy decay E(t) vs. r in simulations (o), with N = 50000 IHS, (a) at low inelasticity 
{a = 0.9,0 = 0.245, r c ~ 70), and (b) at high inelasticity (a = 0.6,0 = 0.05, t c ~ 23), 
compared with theoretical prediction (labeled as 1+2+3, solid line), containing homogeneous 
cooling (1), vorticity diffusion (2) and longitudinal diffusion (3). At the latest simulation times 
(t; in Table 1) both systems are in the nonlinear clustering regime (see fig. 0b). 



Figure shows a similar plot of E{t) in a low density system (0 = 0.05) with a rather large 
inelasticity (a = 0.6), and consequently a short crossover time r c = 23. There is very good 
agreement between theory and simulations over the whole time interval, until the time of the 
latest measurement 77 = 80. At still larger inelasticities and higher densities (a — 0.4, tfi = 0.4) 
the theoretical prediction is an order of magnitude smaller than the simulated E(t). 

Table 1 also shows that the deviations at a = 0.6 between theory and simulations 
increase with increasing density (see <f> — 0.05,0.11,0.4), suggesting that the increase of 
A(to, a)/A(m, 1) with m = (||,-L) as a — * is much larger at high than at low density. A 
possible explanation might be that the Enskog theory for a dense IHS fluid, which is also based 
on molecular chaos, breaks down at higher inelasticities. The effects of dynamic correlations 
are expected to increase strongly with decreasing a, when frequent inelastic collisions force 
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a 





Tc 




Tl 




0.975 


0.245 


340 


27 


500 


0.65* 


0.9 


0.245 


70 


240 


160 


0.8 


0.9 


0.4 


68 


56 


160 


1 


0.6 


0.05 


23 


10 4 


80 


1 


0.6 


0.11 


19 


3000 


160 


1.4 


0.6 


0.4 


15 


150 


80 


4 


0.4 


0.4 


11 


190 


80 


10 



Table I. - Comparison E aim (rt) with E t h eo r{Ti) (*N = 2 x 10 4 ; in all others cases N = 5 x 10 4 ). 



the particles to move more and more parallel. But so far ring kinetic theory |7j for calculating 
transport coefficients in IHS fluids has not yet been developed. 

The present theory is only applicable to thcrmodynamically large systems, because in 
deriving (||) the k-sums have been replaced by k-integrals. In the limit of small inelasticities 
(70 - * 0), the smallest wave number is k ~ 1/L, and the restriction k -C 70 might be violated 
at small L. 

For finite systems deviations from these predictions may occur because of interference effects 
through the periodic boundaries. In elastic hard sphere fluids such effects can be expected for 
times larger than the acoustic traversal time of the system. Here the fastest spreading mode is 
the longitudinal diffusion, and one may expect interference effects when the relevant diffusion 
length, £|| = «/8A||T, satisfies ^ > |L. This implies r > t\\ = L 2 /32A||. For the simulations 
in fig. [j] and Qa the time of the latest measurements, t; = 160, is of the order of T|| = 240, and 
the ratio R(n) = E slm /E theor ~ 0.8. In fig. |b, 77 = 80 < t]| = 10 4 , and ry = 23, and the 
ratio R(ti) = 1. 

In summary, the mode coupling theory, applied to a freely evolving IHS fluid at small 
inelasticity e, explains how the hydrodynamic modes in the elastic range (e <C k < 1) 
are responsible for Haff's homogeneous cooling law (Q) at times t < to/e, and those in the 
dissipative range (k -C e) for the diffusive decay, E(t) ~ r~ d / 2 valid for t < to/e. In fact, 
the full range of hydrodynamic times (t > to) can be evaluated numerically by applying the 
theory of refs. [3a,b] over the full range of hydrodynamic wavenumbers (fc < 1). The solid 
line in figs. ^a,b, representing the sum of homogeneous cooling and transverse and longitudinal 
velocity diffusion, is simply an interpolation between the theoretical long and short time results. 
The agreement between theory and simulations is excellent over the large time interval from 
homogeneous cooling until far into the nonlinear clustering regime. 

We also observe that the motion of the IHS fluid at the largest r- values (see fig. la,b and 
2a, b) is almost 'frozen in', as the kinetic energy has typically decayed by a factor of order 
10 -3 . Therefore, we do not expect a later crossover to truly nonlinear regime, which would 
imply that the decay law, derived here, is asymptotic indeed. At large inelasticities (a ^ 0.4) 
the present theory, which is based on a separation of time scales for small e, is not applicable, 
and the value of E(t), observed in the simulations, is an order of magnitude larger than the 
prediction of the present theory. 
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